--- title: "27_04_23_Model_MH" subtitle: "Modeling forecast count for MH" execute: message: false echo: true warning: false error: true format: html: self-contained: true code-line-numbers: true code-fold: true code-tools: true custom_title_page: false filters: - lightbox lightbox: auto --- ### Log > Version 2023-03-29. Code was initially made in SAS and send by ___Jemisha Apajee___. File Name: "Analysis_common_diseases_step_by_step_20221123.sas" > Version 2023-04-19 Translated to R by ___Javier Silva-Valencia___, adapted to Mental Health paper > Version 2023-05-17 Whitout social problems > Version 2023-06-07 Code debugged by **Percy Soto-Becerra**. The model's problem was fixed. Diagnostic regression assessment was performed. ### Context/Notes - This is part of a mental health paper, where we are analyzing the amount of MH-related visits among several countries (countries part of INTREPID) in a pre and during pandemic period.- The main variables are: Amount of total visits, amount of visits related to MH, 7-9 categories of MH groups, Tipe of visit(in-person, virtual) all this by month from 2018 to 2021 (Except Perú that has data from 2019-2021)- ___This specific code___ is for creating a model that fits the time series and creating a forecast to compare prepandemic and pandemic trend.## Getting ready ### Load packages ```{r, results='hide', warning=FALSE, message = FALSE, collapse = TRUE} if (! require ("pacman" )) install.packages ("pacman" ):: p_load (rio, source (here ("source" , "diag_glmm.R" ))source (here ("source" , "sti_plotter_pre.R" ))source (here ("source" , "sti_model.R" ))source (here ("source" , "sti_model_china.R" ))source (here ("source" , "sti_model_norway.R" ))source (here ("source" , "sti_model_usa.R" ))source (here ("source" , "sti_effect_ploter.R" ))# rm(list = ls()) #To clear the environment and start from zero. # # Colors' palette <- c ("#999999" , # 1 "#E69F00" , # 2 "#56B4E9" , # 3 "#009E73" , # 4 "#F0E442" , # 5 "#0072B2" , # 6 "#D55E00" , # 7 "#CC79A7" ) # 8 ``` ### Read in data ```{r} <- import (here ("Data" , "All_countries_MHcounts.xlsx" ), sheet = "MH_counts" )``` ### View data structure ```{r} datatable (head (mh_visits, 100 )) #Seeing first 10 rows of the dataset ``` ## Cleaning ### Creating variables ```{r} <- mh_visits %>% mutate (date = as.Date (paste (month, "01" , year), format = "%m %d %Y" )) %>% arrange (country, date)``` ```{r} # Select variable of interest <- mh_visits |> :: select (country, #Sorting <- mh_visits1 %>% arrange (country, date, mh_category)# Rename columns <- mh_visits1 %>% rename (numerator = total_counts, #MH counts denominator = den_total) #Total visits per month ``` ### Check for missings ```{r} %>% tabyl (country, period, show_na = TRUE ) #No missings ``` ### Sort categories of MH ```{r} #Eliminando categorias eliminadas <- mh_visits1 %>% filter (! (mh_category %in% c ("Social Problems" ))) # Establecer el nuevo orden para las categorías <- c ("Anxiety and Mood Disorders" , "Bipolar, Schizophrenia and other Psycotic Disorders" ,"Dementia" , "ADHD and Eating Disorders" , "Sleep Disorders" , "Substance-Related and Addictive Disorders" )# Reorder categories <- mh_visits2 %>% filter (mh_category %in% new_order) %>% mutate (mh_category = fct_relevel (factor (mh_category), new_order))``` ### Set up data for modelling ```{r} <- as.Date ("2020-04-01" ) #For all the countries the pandemic start april 2020 *To revise <- as.Date ("2020-02-01" ) #For china pandemic start at feb 2020 <- mh_visits2 %>% # create t = 1, ..., 48 for each country (except peru) mutate (tn = if_else (country != "Peru" , 12 * (year - 2018 )) + month, 12 * (year - 2019 )) + month),time = tn, # for random effects in the residuals #Spliting N# visits into two variables (pre and post start of pandemic) y_count = case_when (!= "China" & date < pandemic_start ~ numerator,== "China" & date < pandemic_start_china ~ numerator,TRUE ~ NA_real_ observed_count = case_when (!= "China" & date >= pandemic_start ~ numerator,== "China" & date >= pandemic_start_china ~ numerator,TRUE ~ NA_real_ # Dummy variables for China (changes in their health system administration) China_Jan19 = if_else (date >= as.Date ("2019-01-01" ), 1 , 0 ),China_Dec20 = if_else (date >= as.Date ("2020-12-01" ), 1 , 0 ),Norway_Feb21 = if_else (time >= 38 , 1 , 0 ), USA_systdown = if_else (time >= 37 & time <= 39 , 1 , 0 ), # Pre-post tag variable country = as.factor (country),pre_post = if_else (period == "pre-pandemic" , 0 , 1 ), pre_postTo = if_else (is.na (y_count), "post" , "pre" ), #pre/post binary variable #"tn" as a factor (will be a classification variable) tn = as.factor (tn), #Creating variable for interaction pre-post exposure (wahsout would be on pre exposure) To = ifelse (period %in% c ("pre-pandemic" , "washout" ), 0 , 1 ), # Creating rate variable for graphics rate = 10000 * numerator / denominator<- mh_visits5 %>% group_by (country, date) %>% summarize (year = max (year),month = max (month),denominator = max (denominator),numerator = sum (numerator),observed_count = sum (observed_count),To = max (To), time = max (time),period = max (period), China_Jan19 = max (China_Jan19), China_Dec20 = max (China_Dec20)%>% ungroup () |> mutate (mh_category = "All common Mental Health disorders" , Norway_Feb21 = if_else (time >= 38 , 1 , 0 ), USA_systdown = if_else (time >= 37 & time <= 39 , 1 , 0 ))|> bind_rows (mh_visits5) -> mh_visits5b``` ### Pre-pandemic behaviour of outcomes (blinded to post-pandemic) ```{r} #| warning: false #| message: false <- min (mh_visits5$ date) #- month(1) <- max (mh_visits5$ date) #+ month(1) <- mh_visits5 |> count (country) |> pull (country) |> as.character ()<- mh_visits5 |> count (mh_category) |> pull (mh_category) |> as.character ()``` ```{r} #| warning: false #| message: false <- expand.grid (country = countryc, mh_category = mh_categoryc)<- combinations %>% mutate (plot = map2 (country, mh_category, ~ sti_plotter_pre (.x, .y, blind = FALSE , data = mh_visits5))) %>% pwalk (~ ggsave (filename = paste0 (..1 , "_" , ..2 , ".png" ), plot = ..3 , device = "png" , path = here ("img" , "exploratory_plots" ), scale = 2 , width = 7 , height = 7 , units = "cm" , dpi = 900 # Extrae la lista de gráficos <- combinations$ plot``` ### Anxiety/Mood ### Argentina ```{r} <- 1 ``` ### Australia ```{r} <- i + 1 ``` ### Canada ```{r} <- i + 1 ``` ### China ```{r} <- i + 1 ``` ### Norway ```{r} <- i + 1 ``` ### Peru ```{r} <- i + 1 ``` ### Singapore ```{r} <- i + 1 ``` ### Sweden ```{r} <- i + 1 ``` ### USA ```{r} <- i + 1 ``` ### Bipolar/Schiz./Psycotic ### Argentina ```{r} <- i + 1 ``` ### Australia ```{r} <- i + 1 ``` ### Canada ```{r} <- i + 1 ``` ### China ```{r} <- i + 1 ``` ### Norway ```{r} <- i + 1 ``` ### Peru ```{r} <- i + 1 ``` ### Singapore ```{r} <- i + 1 ``` ### Sweden ```{r} <- i + 1 ``` ### USA ```{r} <- i + 1 ``` ### Dementia ### Argentina ```{r} <- i + 1 ``` ### Australia ```{r} <- i + 1 ``` ### Canada ```{r} <- i + 1 ``` ### China ```{r} <- i + 1 ``` ### Norway ```{r} <- i + 1 ``` ### Peru ```{r} <- i + 1 ``` ### Singapore ```{r} <- i + 1 ``` ### Sweden ```{r} <- i + 1 ``` ### USA ```{r} <- i + 1 ``` ### ADHD/Eating ### Argentina ```{r} <- i + 1 ``` ### Australia ```{r} <- i + 1 ``` ### Canada ```{r} <- i + 1 ``` ### China ```{r} <- i + 1 ``` ### Norway ```{r} <- i + 1 ``` ### Peru ```{r} <- i + 1 ``` ### Singapore ```{r} <- i + 1 ``` ### Sweden ```{r} <- i + 1 ``` ### USA ```{r} <- i + 1 ``` ### Sleep ### Argentina ```{r} <- i + 1 ``` ### Australia ```{r} <- i + 1 ``` ### Canada ```{r} <- i + 1 ``` ### China ```{r} <- i + 1 ``` ### Norway ```{r} <- i + 1 ``` ### Peru ```{r} <- i + 1 ``` ### Singapore ```{r} <- i + 1 ``` ### Sweden ```{r} <- i + 1 ``` ### USA ```{r} <- i + 1 ``` ### Substance-Related/Addictive ### Argentina ```{r} <- i + 1 ``` ### Australia ```{r} <- i + 1 ``` ### Canada ```{r} <- i + 1 ``` ### China ```{r} <- i + 1 ``` ### Norway ```{r} <- i + 1 ``` ### Peru ```{r} <- i + 1 ``` ### Singapore ```{r} <- i + 1 ``` ### Sweden ```{r} <- i + 1 ``` ### USA ```{r} <- i + 1 ``` ## Creating and Testing the Model ```{r} # Countries and categories <- mh_visits5b |> count (country) |> pull (country) |> as.character ()<- rep (countryc, 7 )<- mh_visits5b |> mutate (mh_category = factor (mh_category, levels = c ("All common Mental Health disorders" , "Anxiety and Mood Disorders" , "ADHD and Eating Disorders" , "Bipolar, Schizophrenia and other Psycotic Disorders" , "Dementia" , "Sleep Disorders" , "Substance-Related and Addictive Disorders" ))) |> count (mh_category) |> pull (mh_category) |> as.character ()``` ## Modelling ### Set global options ```{r} # set options of graphics <- 0.5 <- 0.25 <- c ("#99B933" )<- c ("#002D72" )<- c ("red" )<- c ("skyblue" )<- c ("orange" )<- "2018-01-01" <- "2022-01-01" ``` ### **Outcome:** All common Mental Health disorders ### Countries 1. Argentina2. Australia3. Canada4. China5. Norway6. Peru7. Singapore8. Sweden9. USA### Argentina ```{r} <- 0 ``` ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 5 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot1``` ### Australia ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot2``` ### Canada ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot3``` ### China ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_china (country = countryc[i], mh_category = mh_categoryc[j], p = 2 , q = 1 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot4``` ## Notes - Report in text (maybe in supplementary table?) that including period of shock (t > Dec 2020) does not significantly change the results.### Norway ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_norway (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE , mean_cond = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot5``` ### Peru ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 8 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot6``` ### Singapore ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot7``` ### Sweden ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot8``` ### USA ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` ```{r} # Model <- sti_model_usa (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot9``` #### Table ```{r} # Cargar la biblioteca <- tab |> arrange (country, effect) |> mutate (country_effect = if_else (effect == "Level change" , country, NA_character_ ),effect = case_when (== "Level change" ~ " Level change" ,== "Trend change" ~ " Trend change" |> :: select (- mh_category, - country) |> :: select (country_effect, effect, estimate)# Crear la tabla flextable <- flextable (tab_f)# Cambiar los nombres de las columnas <- set_header_labels (ft, country_effect = "Country" ,effect = "Effect measure" , estimate = "IRR (95%CI), p-value" )# Formatear la tabla <- ft %>% align (align = "left" ) %>% padding (padding = 5 ) %>% fontsize (size = 11 ) %>% bold (part = "header" ) %>% bold (i = which (! is.na (tab_f$ country_effect)), j = "country_effect" ) %>% autofit ()# Guardar como documento de Word save_as_docx (ft, path = here ("table_effects" , paste0 ("table_" , mh_categoryc[j], ".docx" )))``` #### Plot ```{r} #| warning: false #| message: false | plot2 | plot3) / | plot5 | plot6) / | plot8 | plot9)) + plot_layout (guides = 'collect' ) + plot_annotation (tag_levels = 'A' ) & theme (legend.position= "bottom" ) -> plot_effect_finalggsave (filename = paste0 (mh_categoryc[j], ".png" ), plot = plot_effect_final, device = "png" , path = here ("img" , "plot_effects" ), scale = 1.5 , width = 18 , height = 18 , units = "cm" , dpi = 900 )``` ```{r} :: include_graphics (here ("img" , "plot_effects" , paste0 (mh_categoryc[j], ".png" )))``` ### **Outcome:** Anxiety and Mood Disorders ### Countries 1. Argentina2. Australia3. Canada4. China5. Norway6. Peru7. Singapore8. Sweden9. USA### Argentina ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 5 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot1``` ### Australia ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot2``` ### Canada ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot3``` ### China ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_china (country = countryc[i], mh_category = mh_categoryc[j], p = 2 , q = 1 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot4``` ### Norway ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_norway (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE , mean_cond = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot5``` ### Peru ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 8 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot6``` ### Singapore ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot7``` ### Sweden ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot8``` ### USA ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_usa (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot9``` #### Table ```{r} # Cargar la biblioteca <- tab |> arrange (country, effect) |> mutate (country_effect = if_else (effect == "Level change" , country, NA_character_ ),effect = case_when (== "Level change" ~ " Level change" ,== "Trend change" ~ " Trend change" |> :: select (- mh_category, - country) |> :: select (country_effect, effect, estimate)# Crear la tabla flextable <- flextable (tab_f)# Cambiar los nombres de las columnas <- set_header_labels (ft, country_effect = "Country" ,effect = "Effect measure" , estimate = "IRR (95%CI), p-value" )# Formatear la tabla <- ft %>% align (align = "left" ) %>% padding (padding = 5 ) %>% fontsize (size = 11 ) %>% bold (part = "header" ) %>% bold (i = which (! is.na (tab_f$ country_effect)), j = "country_effect" ) %>% autofit ()# Guardar como documento de Word save_as_docx (ft, path = here ("table_effects" , paste0 ("table_" , mh_categoryc[j], ".docx" )))``` #### Plot ```{r} #| warning: false #| message: false | plot2 | plot3) / | plot5 | plot6) / | plot8 | plot9)) + plot_layout (guides = 'collect' ) + plot_annotation (tag_levels = 'A' ) & theme (legend.position= "bottom" ) -> plot_effect_finalggsave (filename = paste0 (mh_categoryc[j], ".png" ), plot = plot_effect_final, device = "png" , path = here ("img" , "plot_effects" ), scale = 1.5 , width = 18 , height = 18 , units = "cm" , dpi = 900 )``` ```{r} :: include_graphics (here ("img" , "plot_effects" , paste0 (mh_categoryc[j], ".png" )))``` ### **Outcome:** ADHD and Eating Disorders ### Countries 1. Argentina2. Australia3. Canada4. China5. Norway6. Peru7. Singapore8. Sweden9. USA### Argentina ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 5 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot1``` ### Australia ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot2``` ### Canada ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot3``` ### China ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_china (country = countryc[i], mh_category = mh_categoryc[j], p = 2 , q = 1 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot4``` ### Norway ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_norway (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE , mean_cond = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot5``` ### Peru ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 8 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot6``` ### Singapore ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` `r mh_categoryc[j]` for Singapore is not appropiate due to missing data in numerator.```{r} <- data.frame (x = seq (1 , 48 , 1 ), y = seq (1 , 48 , 1 ))|> ggplot (aes (x = x, y = y)) + theme_void () + annotate (geom = "text" , x = 24 , y = 24 , label = "Not calculable" , color = "black" , size = 6 ) -> plot7``` ### Sweden ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot8``` ### USA ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_usa (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot9``` #### Table ```{r} # Cargar la biblioteca <- tab |> arrange (country, effect) |> mutate (country_effect = if_else (effect == "Level change" , country, NA_character_ ),effect = case_when (== "Level change" ~ " Level change" ,== "Trend change" ~ " Trend change" |> :: select (- mh_category, - country) |> :: select (country_effect, effect, estimate)# Crear la tabla flextable <- flextable (tab_f)# Cambiar los nombres de las columnas <- set_header_labels (ft, country_effect = "Country" ,effect = "Effect measure" , estimate = "IRR (95%CI), p-value" )# Formatear la tabla <- ft %>% align (align = "left" ) %>% padding (padding = 5 ) %>% fontsize (size = 11 ) %>% bold (part = "header" ) %>% bold (i = which (! is.na (tab_f$ country_effect)), j = "country_effect" ) %>% autofit ()# Guardar como documento de Word save_as_docx (ft, path = here ("table_effects" , paste0 ("table_" , mh_categoryc[j], ".docx" )))``` #### Plot ```{r} #| warning: false #| message: false | plot2 | plot3) / | plot5 | plot6) / | plot8 | plot9)) + plot_layout (guides = 'collect' ) + plot_annotation (tag_levels = 'A' ) & theme (legend.position= "bottom" ) -> plot_effect_finalggsave (filename = paste0 (mh_categoryc[j], ".png" ), plot = plot_effect_final, device = "png" , path = here ("img" , "plot_effects" ), scale = 1.5 , width = 18 , height = 18 , units = "cm" , dpi = 900 )``` ```{r} :: include_graphics (here ("img" , "plot_effects" , paste0 (mh_categoryc[j], ".png" )))``` ### **Outcome:** Bipolar, Schizophrenia and other Psycotic Disorders ### Countries 1. Argentina2. Australia3. Canada4. China5. Norway6. Peru7. Singapore8. Sweden9. USA### Argentina ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 5 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot1``` ### Australia ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot2``` ### Canada ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot3``` ### China ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_china (country = countryc[i], mh_category = mh_categoryc[j], p = 2 , q = 1 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot4``` ### Norway ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_norway (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE , mean_cond = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot5``` ### Peru ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 8 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot6``` ### Singapore ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot7``` ### Sweden ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot8``` ### USA ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_usa (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot9``` #### Table ```{r} # Cargar la biblioteca <- tab |> arrange (country, effect) |> mutate (country_effect = if_else (effect == "Level change" , country, NA_character_ ),effect = case_when (== "Level change" ~ " Level change" ,== "Trend change" ~ " Trend change" |> :: select (- mh_category, - country) |> :: select (country_effect, effect, estimate)# Crear la tabla flextable <- flextable (tab_f)# Cambiar los nombres de las columnas <- set_header_labels (ft, country_effect = "Country" ,effect = "Effect measure" , estimate = "IRR (95%CI), p-value" )# Formatear la tabla <- ft %>% align (align = "left" ) %>% padding (padding = 5 ) %>% fontsize (size = 11 ) %>% bold (part = "header" ) %>% bold (i = which (! is.na (tab_f$ country_effect)), j = "country_effect" ) %>% autofit ()# Guardar como documento de Word save_as_docx (ft, path = here ("table_effects" , paste0 ("table_" , mh_categoryc[j], ".docx" )))``` #### Plot ```{r} #| warning: false #| message: false | plot2 | plot3) / | plot5 | plot6) / | plot8 | plot9)) + plot_layout (guides = 'collect' ) + plot_annotation (tag_levels = 'A' ) & theme (legend.position= "bottom" ) -> plot_effect_finalggsave (filename = paste0 (mh_categoryc[j], ".png" ), plot = plot_effect_final, device = "png" , path = here ("img" , "plot_effects" ), scale = 1.5 , width = 18 , height = 18 , units = "cm" , dpi = 900 )``` ```{r} :: include_graphics (here ("img" , "plot_effects" , paste0 (mh_categoryc[j], ".png" )))``` ### **Outcome:** Dementia ### Countries 1. Argentina2. Australia3. Canada4. China5. Norway6. Peru7. Singapore8. Sweden9. USA### Argentina ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 5 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot1``` ### Australia ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot2``` ### Canada ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot3``` ### China ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_china (country = countryc[i], mh_category = mh_categoryc[j], p = 2 , q = 1 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot4``` ### Norway ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_norway (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE , mean_cond = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot5``` ### Peru ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 8 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot6``` ### Singapore ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot7``` ### Sweden ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot8``` ### USA ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_usa (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot9``` #### Table ```{r} # Cargar la biblioteca <- tab |> arrange (country, effect) |> mutate (country_effect = if_else (effect == "Level change" , country, NA_character_ ),effect = case_when (== "Level change" ~ " Level change" ,== "Trend change" ~ " Trend change" |> :: select (- mh_category, - country) |> :: select (country_effect, effect, estimate)# Crear la tabla flextable <- flextable (tab_f)# Cambiar los nombres de las columnas <- set_header_labels (ft, country_effect = "Country" ,effect = "Effect measure" , estimate = "IRR (95%CI), p-value" )# Formatear la tabla <- ft %>% align (align = "left" ) %>% padding (padding = 5 ) %>% fontsize (size = 11 ) %>% bold (part = "header" ) %>% bold (i = which (! is.na (tab_f$ country_effect)), j = "country_effect" ) %>% autofit ()# Guardar como documento de Word save_as_docx (ft, path = here ("table_effects" , paste0 ("table_" , mh_categoryc[j], ".docx" )))``` #### Plot ```{r} #| warning: false #| message: false | plot2 | plot3) / | plot5 | plot6) / | plot8 | plot9)) + plot_layout (guides = 'collect' ) + plot_annotation (tag_levels = 'A' ) & theme (legend.position= "bottom" ) -> plot_effect_finalggsave (filename = paste0 (mh_categoryc[j], ".png" ), plot = plot_effect_final, device = "png" , path = here ("img" , "plot_effects" ), scale = 1.5 , width = 18 , height = 18 , units = "cm" , dpi = 900 )``` ```{r} :: include_graphics (here ("img" , "plot_effects" , paste0 (mh_categoryc[j], ".png" )))``` ### **Outcome:** Sleep Disorders ### Countries 1. Argentina2. Australia3. Canada4. China5. Norway6. Peru7. Singapore8. Sweden9. USA### Argentina ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 5 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot1``` ### Australia ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot2``` ### Canada ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot3``` ### China ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_china (country = countryc[i], mh_category = mh_categoryc[j], p = 2 , q = 1 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot4``` ### Norway ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_norway (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE , mean_cond = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot5``` ### Peru ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 8 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot6``` ### Singapore ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` `r mh_categoryc[j]` for Singapore is not appropiate due to missing data in numerator.```{r} <- data.frame (x = seq (1 , 48 , 1 ), y = seq (1 , 48 , 1 ))|> ggplot (aes (x = x, y = y)) + theme_void () + annotate (geom = "text" , x = 24 , y = 24 , label = "Not calculable" , color = "black" , size = 6 ) -> plot7``` ### Sweden ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot8``` ### USA ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_usa (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot9``` #### Table ```{r} # Cargar la biblioteca <- tab |> arrange (country, effect) |> mutate (country_effect = if_else (effect == "Level change" , country, NA_character_ ),effect = case_when (== "Level change" ~ " Level change" ,== "Trend change" ~ " Trend change" |> :: select (- mh_category, - country) |> :: select (country_effect, effect, estimate)# Crear la tabla flextable <- flextable (tab_f)# Cambiar los nombres de las columnas <- set_header_labels (ft, country_effect = "Country" ,effect = "Effect measure" , estimate = "IRR (95%CI), p-value" )# Formatear la tabla <- ft %>% align (align = "left" ) %>% padding (padding = 5 ) %>% fontsize (size = 11 ) %>% bold (part = "header" ) %>% bold (i = which (! is.na (tab_f$ country_effect)), j = "country_effect" ) %>% autofit ()# Guardar como documento de Word save_as_docx (ft, path = here ("table_effects" , paste0 ("table_" , mh_categoryc[j], ".docx" )))``` #### Plot ```{r} #| warning: false #| message: false | plot2 | plot3) / | plot5 | plot6) / | plot8 | plot9)) + plot_layout (guides = 'collect' ) + plot_annotation (tag_levels = 'A' ) & theme (legend.position= "bottom" ) -> plot_effect_finalggsave (filename = paste0 (mh_categoryc[j], ".png" ), plot = plot_effect_final, device = "png" , path = here ("img" , "plot_effects" ), scale = 1.5 , width = 18 , height = 18 , units = "cm" , dpi = 900 )``` ```{r} :: include_graphics (here ("img" , "plot_effects" , paste0 (mh_categoryc[j], ".png" )))``` ### **Outcome:** Substance-Related and Addictive Disorders ### Countries 1. Argentina2. Australia3. Canada4. China5. Norway6. Peru7. Singapore8. Sweden9. USA### Argentina ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 5 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot1``` ### Australia ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot2``` ### Canada ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot3``` ### China ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_china (country = countryc[i], mh_category = mh_categoryc[j], p = 2 , q = 1 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot4``` ### Norway ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_norway (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE , mean_cond = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot5``` ### Peru ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 8 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot6``` ### Singapore ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 11 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot7``` ### Sweden ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model (country = countryc[i], mh_category = mh_categoryc[j], p = 6 , q = 0 , data = mh_visits5b)``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot8``` ### USA ```{r} #| message: false #| warning: false # Set the counter <- i + 1 if (i %% 9 == 0 ) {<- (i - 1 ) %/% 9 + 1 else {<- i %/% 9 + 1 ``` - **Country:** `r countryc[i]` - **Mental Health Category:** `r mh_categoryc[j]` ```{r} # Model <- sti_model_usa (country = countryc[i], mh_category = mh_categoryc[j], p = 1 , q = 0 , data = mh_visits5b, shock_exclude = TRUE )``` - The autocorrelation function of residuals from a naive glm that does not account for autocorrelation is the following: ```{r} #| message: false #| warning: false $ auto_plot; mod_full$ arima``` - The normalized residuasl from a GLMMM that take into account dependence of residuls around the time shows no autocorrelation:```{r} #| message: false #| warning: false $ dx.plot``` ```{r} summary (mod_full$ mod)``` ```{r} :: model_parameters (mod_full$ mod) |> :: select (Parameter, Coefficient, CI_low, CI_high, p) |> filter (Parameter %in% c ("To" , "To:I(time - tpand)" )) |> as_data_frame () |> rename (est = Coefficient, ll = CI_low, ul = CI_high) |> mutate (effect = c ("Level change" , "Trend change" ), across (c (est, ll, ul), ~ (round (exp (.x) , 2 ))), p = if_else (p < 0.001 , "p < 0.001" , as.character (paste0 ("p = " , round (p, 3 )))), estimate = str_glue ("{est} ({ll} to {ul}), {p}" ), estimate = as.character (estimate), mh_category = mh_categoryc[j], country = countryc[i]|> :: select (country, mh_category, effect, estimate) -> tab0<- tab0 |> bind_rows (tab)|> :: kable ()``` ```{r} #| warning: false #| message: false sti_effect_ploter (mod_full, alfa = alfa0,alfa_ribbon = alfa_ribbon0,color_point = color_point0, color_line = color_line0, color_ribbon_ex = color_ribbon_ex0,color_ribbon_def = color_ribbon_def0, color_interrupt = color_interrupt0,mindate = mindate0, maxdate = maxdate0) -> plot9``` #### Table ```{r} # Cargar la biblioteca <- tab |> arrange (country, effect) |> mutate (country_effect = if_else (effect == "Level change" , country, NA_character_ ),effect = case_when (== "Level change" ~ " Level change" ,== "Trend change" ~ " Trend change" |> :: select (- mh_category, - country) |> :: select (country_effect, effect, estimate)# Crear la tabla flextable <- flextable (tab_f)# Cambiar los nombres de las columnas <- set_header_labels (ft, country_effect = "Country" ,effect = "Effect measure" , estimate = "IRR (95%CI), p-value" )# Formatear la tabla <- ft %>% align (align = "left" ) %>% padding (padding = 5 ) %>% fontsize (size = 11 ) %>% bold (part = "header" ) %>% bold (i = which (! is.na (tab_f$ country_effect)), j = "country_effect" ) %>% autofit ()# Guardar como documento de Word save_as_docx (ft, path = here ("table_effects" , paste0 ("table_" , mh_categoryc[j], ".docx" )))``` #### Plot ```{r} #| warning: false #| message: false | plot2 | plot3) / | plot5 | plot6) / | plot8 | plot9)) + plot_layout (guides = 'collect' ) + plot_annotation (tag_levels = 'A' ) & theme (legend.position= "bottom" ) -> plot_effect_finalggsave (filename = paste0 (mh_categoryc[j], ".png" ), plot = plot_effect_final, device = "png" , path = here ("img" , "plot_effects" ), scale = 1.5 , width = 18 , height = 18 , units = "cm" , dpi = 900 )``` ```{r} :: include_graphics (here ("img" , "plot_effects" , paste0 (mh_categoryc[j], ".png" )))``` ## Reproducibility Ticket ```{r} sessionInfo ()``` 
